load "/data/data_s11_3/jlli/ARPS/Entrainment/resolvable/functions/unpack_16bit.ncl"

external QVFLUX "/data/data_s11_3/jlli/ARPS/Entrainment/resolvable_new/fortran/qvflux_interp99.so"

begin

; Generate the ARPS data file names
second_tmp  = ispan(1210, 1980, 10)
second_str  = tostring(second_tmp)
second      = new(dimsizes(second_tmp), string)
num_second  = dimsizes(second_tmp)
second(:)= "00"+second_str(:)

arps_file     = "newfric_50m_cd1-2A6.hdf"+second
dir_arpsdata  = "/data/data_s11_3/jlli/ARPS/LiJinlin/new_friction_runs_100m_WK1982_dewp8wsm1wbc_ns_cdm2_61_totalvg/newfric_100m_cd1-2A6.outdata-10second/newfric_50m_cd1-2A6.data"+second+"/"
dir_arpsbas   = "/data/data_s11_3/jlli/ARPS/LiJinlin/new_friction_runs_100m_WK1982_dewp8wsm1wbc_ns_cdm2_61_totalvg/newfric_100m_cd1-2A6.outdata-10second/newfric_50m_cd1-2A6.data000000/newfric_50m_cd1-2A6.hdfgrdbas"

; Calculate the altitude of 'Scalar Value'
delim      = " "
datatsall  = readAsciiTable("/data/data_s11_3/jlli/ARPS/Entrainment/sounding/CN_61cases_Vg_sounding_APRS_modelleve_nomove_ns.txt", 1, "string", 2)
datats     = datatsall(0:81,0)
level_tmp  = tofloat(str_get_field(datats, 1, delim))
levels     = level_tmp(::-1)

level = new(82, float)
do iz = 0, 80
    level(iz+1) = (levels(iz)+levels(iz+1))/2.0
end do
level(0) = -20.0


; Parameter Definition
nsub = 9

nz = 83
ny = 903
nx = 903

skipy_s = 300
skipy_e = 550
skipx_s = 300
skipx_e = 550
eallys = skipy_e-skipy_s
eallxs = skipx_e-skipx_s
; eallz = 49 ; 5km
eallzs = 56 ; 7km
; eallz = 78 ; 15km

eallz = eallzs + 1
eally = eallys + 1
eallx = eallxs + 1

eallz_news = eallzs*nsub
eally_news = eallys*nsub
eallx_news = eallxs*nsub

eallz_new = eallz_news + 1
eally_new = eally_news + 1
eallx_new = eallx_news + 1

dx = 100.0
dy = 100.0
wtrhold = 2.0
qtrhold = -1.0 ; it means that the influence of cloud is not considerd when this value is less than zero


rho = new(eallzs, float)
qvinflux_x_iz = new((/num_second,eallzs/), float)
qvinflux_y_iz = new((/num_second,eallzs/), float)
qvinflux_z_iz = new((/num_second,eallzs/), float)
qvinflux_iz = new((/num_second,eallzs/), float)
qvoutflux_x_iz = new((/num_second,eallzs/), float)
qvoutflux_y_iz = new((/num_second,eallzs/), float)
qvoutflux_z_iz = new((/num_second,eallzs/), float)
qvoutflux_iz = new((/num_second,eallzs/), float)
qvflux_e_iz = new((/num_second,eallzs/), float)
qvflux_d_iz = new((/num_second,eallzs/), float)
v_iz = new((/num_second,eallzs/), float)

; Calculate entrainment rate =========================================================================================================

; Read rhobar data
arpsbasdata  = addfile(dir_arpsbas, "r")
rhobar_tmp   = arpsbasdata->rhobar  ; Base state air density
rhobar       = unpack_16bit(rhobar_tmp)
rho(:)       = rhobar(0:eallzs-1,449,449)
level_sub    = level(0:eallz-1)

do it = 0, num_second-1
    print("The current time is '"+second_tmp(it)+"'")
    arpsdata  = addfile(dir_arpsdata(it)+arps_file(it), "r")
    uwind_tmp = arpsdata->u  ; U component of total velocity
    vwind_tmp = arpsdata->v  ; V component of total velocity
    wwind_tmp = arpsdata->w  ; Vertical velocity
    qv_tmp    = arpsdata->qv ; Water vapor specific humidity
    qc_tmp    = arpsdata->qc ; Cloud water
    qi_tmp    = arpsdata->qi ; Cloud ice
    ; qr_tmp    = arpsdata->qr ; Rain water
    ; qs_tmp    = arpsdata->qs ; Snow
    ; qh_tmp    = arpsdata->qh ; Hail

    uwind     = unpack_16bit(uwind_tmp)
    vwind     = unpack_16bit(vwind_tmp)
    wwind     = unpack_16bit(wwind_tmp)
    qv        = unpack_16bit(qv_tmp)
    qc        = unpack_16bit(qc_tmp)
    qi        = unpack_16bit(qi_tmp)
    ; qr        = unpack_16bit(qr_tmp)
    ; qs        = unpack_16bit(qs_tmp)
    ; qh        = unpack_16bit(qh_tmp)

    qc        = qc*1000
    qi        = qi*1000
    ; qr        = qr*1000
    ; qs        = qs*1000
    ; qh        = qh*1000

    tq       = qc+qi ;+qr+qs+qh  ;Total q

    tq_sub  = tq(0:eallzs-1,skipy_s:skipy_e-1,skipx_s:skipx_e-1) ;Total q
    qv_sub  = qv(0:eallzs,skipy_s:skipy_e,skipx_s:skipx_e)
    u_sub   = uwind(0:eallzs-1,skipy_s:skipy_e-1,skipx_s:skipx_e)
    v_sub   = vwind(0:eallzs-1,skipy_s:skipy_e,skipx_s:skipx_e-1)
    w_sub   = wwind(0:eallzs,skipy_s:skipy_e-1,skipx_s:skipx_e-1)

    QVFLUX::qvflux_interp99(nsub,eallz,eally,eallx,eallzs,eallys,eallxs,eallz_new,eally_new,eallx_new,eallz_news,eally_news,eallx_news,\
                            u_sub,v_sub,w_sub,qv_sub,tq_sub,rho,level_sub,dx,dy,wtrhold,qtrhold,\
                            qvinflux_x_iz(it,:),qvinflux_y_iz(it,:),qvinflux_z_iz(it,:),qvinflux_iz(it,:),\
                            qvoutflux_x_iz(it,:),qvoutflux_y_iz(it,:),qvoutflux_z_iz(it,:),qvoutflux_iz(it,:),\
                            qvflux_e_iz(it,:),qvflux_d_iz(it,:),v_iz(it,:))

end do; time

dir_generate = "/data/data_s11_3/jlli/ARPS/Entrainment/resolvable_new/entrainment_qvflux/"

setfileoption("nc", "Format", "NETCDF4")
outfile = addfile(dir_generate+"qvflux_interpolate_100m_dewp8wsm1wbc_totalvg_1210-1980_test_f99_wider1_int9.nc", "c")
outfile->qvinflux_x_iz     = qvinflux_x_iz
outfile->qvinflux_y_iz     = qvinflux_y_iz
outfile->qvinflux_z_iz     = qvinflux_z_iz
outfile->qvinflux_iz       = qvinflux_iz
outfile->qvoutflux_x_iz     = qvoutflux_x_iz
outfile->qvoutflux_y_iz     = qvoutflux_y_iz
outfile->qvoutflux_z_iz     = qvoutflux_z_iz
outfile->qvoutflux_iz       = qvoutflux_iz
outfile->qvflux_e_iz        = qvflux_e_iz
outfile->qvflux_d_iz        = qvflux_d_iz
outfile->v_iz               = v_iz

end